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Abstract 

Previous Monte Carlo investigations by Wojciechowski et al. have found two unusual phases 
in two-dimensional systems of anisotropic hard particles: a tetratic phase of four-fold symmetry 
for hard squares [Comp. Methods in Science and Tech., 10: 235-255, 2004], an d a nonperiodic 
degenerate solid phase for hard-disk dimers [Phys. Rev. Lett., 66: 3168-3171, 1991]. In this work, 
we study a system of hard rectangles of aspect ratio two, i.e., hard-square dimers (or dominos), and 
demonstrate that it exhibits a solid phase with both of these unusual properties. The solid shows 
tetratic, but not nematic, order, and it is nonperiodic having the structure of a random tiling of the 
square lattice with dominos. We obtain similar results with both a classical Monte Carlo method 
using true rectangles and a novel molecular dynamics algorithm employing rectangles with rounded 
corners. It is remarkable that such simple convex two-dimensional shapes can produce such rich 
phase behavior. Although we have not performed exact free-energy calculations, we expect that 
the random domino tiling is thermodynamically stabilized by its degeneracy entropy, well-known 
to be 1.79/cg per particle from previous studies of the dimer problem on the square lattice. Our 
observations are consistent with a KTHNY two-stage phase transition scenario with two continuous 
phase transitions, the first from isotropic to tetratic liquid, and the second from tetratic liquid to 
solid. 



I. INTRODUCTION 



Hard-particle systems have provided a simple and rich model for investigating phase 
behavior and transport in atomic and molecular materials. It is long-known that a pure 
hard-core exclusion potential can lead to a variety of behaviors depending on the degree of 
anisotropy of the particles, including the occurrence of isotropic and nematic liquids, layered 
smectic, and ordered solid phases. 1 Through computer investigations of various particle 
shapes, other phases have been found, such as the biaxial 2 (recently synthesized in the 
laboratory- ) and cubatic phases in three dimensions, in which the axes of symmetry of 
the individual particles align along two or three perpendicular axes (directors). One only 
need look at simple shapes in two dimensions to discover interesting phases. In recent work, 
Wojciechowski et al. studied hard squares and found the first example of a tetratic liquid 
phase at intermediate densities.™ In a tetratic liquid, there is (quasi) long-range orientational 
ordering along two perpendicular axes, but only short-range translational ordering. The solid 
phase is the expected square lattice, with quasi-long-range periodic ordering. On the other 
hand, by studying hard-disk dimers (two disks fused at a point on their boundary), they 
have identified the first example of a nonperiodic solid phase at high densities.— In this 
phase, the centroids of the particles are ordered on the sites of a triangular lattice. However, 
the orientations of the dimers are disordered, leading to a high degeneracy entropy of the 
nonperiodic solid and a lower free energy as compared to periodic solids. 

In this paper, we look at systems of rectangles of aspect ratio a = a/b = 2, i.e., hard- 
square dimers (or dominos). Since the aspect ratio is far from unity, it is not clear a priori 
whether nematic or tetratic orientational ordering (or both) will appear. Recent Density 
Functional Theory calculations,— extending previous work based on scaled-particle theory/ 
have predicted that for a = 2 the tetratic phase is only metastable with respect to the 
ordered solid phase in which all particles are aligned. However, these calculations are only 
approximate and the authors point out that tetratic order is still possible in spatially ordered 
phases. An obvious candidate for forming a stable tetratic phase are dominos: two dominos 
paired along their long edges form a square, and these squares can then form a square lattice 
assuming one of two random orientations, thus forming a tetratic phase with degeneracy 
entropy of ln(\/2). In fact, one does not need to pair up the rectangles but rather simply tile 
a square lattice with dominos which randomly assume one of the two preferred perpendicular 
directions. The degeneracy entropy of this domino tiling has been calculated exactly to 
be (2G/rr)k B w 0.58313A; B , 8 ' 9 where G = J2n=o (-l) n (2n + l)" 2 w 0.91597 is Catalan's 
constant. At high densities, free- volume theory^ predicts that the configurational entropy 
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(per particle) diverges like 

S FY ~ / ln(l - <p/(j) c ) + S con{ , 

where / is the (effective) number of degrees of freedom per particle, C is the volume fraction 
(density) at close packing, and S^onf is an additive constant due to collective exclusion- volume 
effects. Therefore, the densest solid is thermodynamically favored, but if several solids have 
the same density the additive factor matters. Therefore, for hard rectangles, for which the 
maximal density is <ft c = 1 and is achieved by a variety of packings, the degeneracy entropy can 
dominate S^onf and thus the nonperiodic random tiling can be thermodynamically favored. 
Indeed, we demonstrate that our simulations of the hard-domino system produce high-density 
phases with structures very similar to that of a random covering of the square lattice with 
dimers. 

The phase transitions in two-dimensional systems are of interest to the search for contin- 
uous KTHN Y 10 i n i 12 transitions between the disordered liquid and the ordered solid phase. 
At present there is no agreement on the nature of the transition even for the hard-disk sys- 
tem. A previous study of the melting of a square-lattice crystal, stabilized by the addition of 
three-body interactions, found evidence of a (direct) first-order melting.^ Our observations 
for the domino system are relatively consistent with a KTHNY-like two-stage transition: a 
continuous phase transition from an isotropic to a tetratic liquid with (quasi-) long-range 
tetratic order around ~ 0.7, and then another continuous transition from tetratic liquid to 
tetratic solid with quasi-long-range translational order <fi « 0.8. However, we cannot rule out 
the possibility of a weak first-order phase transition between the two phases without more 
detailed simulations. 

This paper is organized as follows. In Section |TlJ we present the simulation techniques 
used to generate equilibrated systems at various densities. In Section IIII1 we analyze the 
properties of the various states, focusing on the orientational and translational ordering in 
the high-density phases. We conclude with a summary of the results and suggestions for 
future work in Section HV1 

II. SIMULATION TECHNIQUES 

In this section, we provide additional details on the MC and MD algorithms we imple- 
mented. It is important to point out that it is essential to implement techniques for speeding 
up the near-neighbour search, in both MC and MD. For rectangles with a small aspect ratio, 
we employ the well-known technique of splitting the domain of simulation into cells (bins) 
larger than a particle diameter D = \J a 2 + b 2 , and consider as neighbors only particles whose 
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centroids belong to neighboring cells. Additional special techniques more suitable for very 
aspherical particles or systems near jamming are described in Ref. 



A. Monte Carlo 

We have implemented a standard MC algorithm in the NVT ensemble, with the additional 
provision of changing the density by growing or shrinking the particles in small increments. 
Each rectangle is described by the location of its centroid (x, y) and orientation 9. For 
increased computational speed the pair (sin 9, cos 9) may be used to represent the orientation. 
In a trial MC step, a rectangle is chosen at random and its coordinates are changed slightly, 
either translationally (Ax, Ay) or orientationally (Ad). Every move has an equal chance of 
being translational or orient at ional. The rectangle's new position is then compared against 
nearby rectangles for overlap; if there is no overlap, the trial move is accepted. We call 
a sequence of N trials a cycle. The simulation evolves through stages, defined by a speed 
n cydes/ stage- At the end of a cycle, pressure data are collected by the virtual-scaling method 
of Eppenga and Frenkel.— Namely, p = PV/NkT = 1 + (pa/2 where a is the rate at which 
growing the particles causes overlaps. At the end of a stage, order parameters and other 
statistics are collected, and then the packing fraction is changed by a small value A0; it 
may be increased, decreased, or not changed at all. If A<p > 0, then <p cannot necessarily 
change by A(f> every stage, because the increase could create overlaps. We scale down the 
increase by factors of 2 until a A0 e g- is found that does not cause any overlaps when applied. 
Typical values for runs are n cyc i es / stage = 1000 and Acf) = ±1 x 10~ 5 . Since there is a limit on 
how fast one can increase the density in such a Monte Carlo simulation, especially at very 
high densities, we use Molecular Dynamics to compress systems to close packing. 

The overlap test is by far the largest computational bottleneck in the MC program. The 
overlap test for two rectangles is based on the following fact: Two rectangles R\ and R2 do 
not overlap if and only if a separating line I can be drawn such that all four corners of R\ 
lie on one side of the line and all four corners of R2 lie on the other side.— The corners of 
both rectangles are allowed to coincide with i. Without loss of generality, we may assume £ 
is drawn parallel to one of the rectangles' major axes and runs exactly along that rectangle's 
side. The problem of testing all possible lines £ thus reduces to testing the eight lines that 
coincide with the edges of R\ and i?2- The test can be optimized somewhat further, as 
illustrated in Fig. ^ An axis a of R\ is chosen. The distance from the center of R2 to a is 
found. Then the distance do of closest approach of R2 to a is found by subtracting a sine and 
a cosine; this distance corresponds to the corner of R 2 that is closest to a. By comparing 
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Figure 1: Illustration of the optimized overlap test for two rectangles. The axes are a and b, with 
semiaxes a and b, and the length from a to the closest corner of R2 is do- (Left) If do > b, then the 
rectangles do not overlap. (Right) If do < b, then the rectangles overlap. 

do with the length b of the other semiaxis of R\, two possible lines £, corresponding to two 
opposite sides of Ri, can be tested at once. If do < b, there is an overlap. In this way, four 
different values of do are calculated; one for each axis of each rectangle. If no comparison 
finds an overlap, there is no overlap. 



B. Molecular Dynamics 



MC simulations are typically the most efficient when one is only interested in stable equi- 
librium properties. We have previously developed a molecular dynamics (MD) algorithm 
aimed at studying hard nonspherical particles and applied it to systems of hard ellipses and 
ellipsoids. 14 We have since generalized the implementation to also handle "superellipses" and 
"superellipsoids", which are generalized smooth convex shapes capable of approximating cen- 
trally symmetric shapes with sharp corners such as rectangles. A superellipse with semiaxes 
a and b is given by the equation 



2C 



+ 



2C 



i/C 



< 1, 



where ( > 1 is an exponent. We add an exponent l/( above in order to properly normalize 
the convex function defining the particle shape, even though it is not strictly necessary. 
When ( = 1 we get the simple ellipse, and when £ — > 00 we obtain a rectangle with sides 
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Figure 2: An snapshot of a few superellipses (exponent £ = 7.5) used in the MD simulations. It 
can be seen that the particle shape is very close to a rectangle. 

2a and 2b. The higher the exponent the sharper the corners become. The smoothing of the 
corners of the rectangle enables us to apply our collision-driven MD algorithm,— with few 
changes from the case of ellipses. Details of this implementation will be given elsewhere. The 
floating-point cost of the algorithm increases as the exponent increases, while the numerical 
stability decreases. We have used an exponent ( = 7.5 for the studies presented here (for this 
exponent the ratio of the areas of the superellipse and the true rectangle is 0.9934). Figure 
12] gives an illustration of the particle shape. 

There are some advantages of the MD simulation over MC. The shapes of the particles can 
change arbitrarily fast in an easily controlled manner by simply adding a dynamic growth 
rate 7 = da/dt = adb/dt. If 7 > 0, i.e., the density is increasing, two colliding particles 
simply get an extra repulsive boost that ensures no overlaps are created. The velocities 
are periodically rescaled to T = 1 to compensate for the induced heating or cooling due 
to the particle growth.— In general, (common) MC methods do not work well near close 
packing, while MD methods, especially event-driven ones, can successfully be used to study 
the neighbourhood of jamming points. Additionally, pressure measurement is more natural 
in the MD method, as the pressure can be directly obtained from time averages of the 
momentum exchange in binary collisions between particles. We have found this pressure 
measurement to be much more precise than using virtual particle scaling in MC simulations. 
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III. RESULTS 



By using either the MC or the MD algorithm with small particle growth rate (A0 or 
7), we have traced the (quasi)equilibrium phase behavior of systems of dominos over a 
range of densities. In this section, we present several techniques for measuring orientational 
and translational order for a given configuration of particles, as well as the results of such 
measurements for the generated states. We have tested our codes by first applying them 
to hard squares and comparing the results to those in Ref. U, and we have observed good 
quantitative agreement throughout. Our MC pressure measurement systematically slightly 
underestimates the pressure compared to the NPT ensemble used in Ref. U and to our MD 
simulations. We present some of the results for the MC, and others for the MD simulations, 
marking any quantitative differences. The two techniques always produced qualitatively 
identical results. 

Describing the statistical properties of the observed states would require specifying all 
of the n-particle correlation functions. The most important is the pair correlation function 
t/ 2 (r, ip, AO). Given a particle, g2{r, ip, AO) is the probability density of finding another parti- 
cle whose centroid is a distance r away (from the centroid of the particle), at a displacement 
angle of ip (relative to the first particle's coordinate axes), and with a orientation of AO 
(relative to the particle's orientation). The normalization of g<i is such that it is identically 
unity for an ideal gas. We will use an equivalent representation where we fix a particle at 
the origin such that the longer rectangle axis is along the x axes, and represent pair cor- 
relations with g2(Ax,Ay,A6), giving the probability density that there is another particle 
whose centroid is at position (Ax, Ay) and whose major axis is at a relative angle of AO. 
Since a three-dimensional function is rather difficult to calculate accurately and visualize, we 
can separate the translational and orientational components and average over some of the 
dimensions to reduce it to a one- or two-dimensional function. 

A. Equation Of State 

The pressure as a function of density can be most accurately measured in the MD simu- 
lations. There is no exact theory that can predict the entire equation of state (EOS) for a 
given many-particle system. However, there are two simple theories that produce remarkably 
good predictions for a variety of systems studied in the literature. For the isotropic fluid 
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(gas) phase of a system of hard dominos, scaled-particle theory (SPT)iI predicts 

P=^- + -^^ (1) 



and modifications to account for possible orientational ordering are discussed in Refs. 
For the solid phase, the free-volume (FV) theory predicts a divergence of the pressure near 
close packing of the form 

P=Y7^7Ti ( 2 ) 



and (liquid-state) density functional theory can be used to make quantitative predictions 
at intermediate densities.— For superellipses with exponent ( = 7.5 the maximal density is 
somewhat less than 1 and we take it to be equal to the ratio of the areas of the particle and 
a true rectangle, C ~ 0.9934. 

The numerical EOS from the NVT MD simulation are shown in Fig. El for both a slow 
compression starting from an isotropic liquid and a decompression starting from a perfect 
random domino tiling generated with the help of random spanning trees, as explained in 
Ref. LL8J. We note that the random domino tiling used was generated inside a square box 
(see Fig. ITT]) even though periodic boundary conditions were used in the actual simulation. 
We expect this to have a very small effect.— It is clearly seen from the figure that there is a 
transition from the liquid to the solid branch in the region m 0.7 and <p ~ 0.8, although no 
clear discontinuities or a hysteresis loop are seen (which would be indicative of a first-order 
phase transition). Compressing an isotropic liquid invariably freezes some defects and thus 
the jamming density is smaller (and the pressure is thus higher) than in the perfect crystal. 



B. Orientational Order 



Orientational order can be measured via the orientational correlation function of order m 

G m (r) = (cos(mA#)) r , (3) 

where m is an integer and the average is taken over all pairs of particles that are at a 
distance between r and r + dr apart from each other. The one-dimensional function G m (r) 
can be thought of as giving normalized Fourier components of the distribution of relative 
orientations versus interparticle distance. When m = 2, it measures the degree of nematic 
ordering (parallel alignment of the particles' major axes), and when m — 4 it measures the 
degree of tetratic ordering (parallel alignment of the particles' axes). The infinite-distance 
value \\mG m {r) = S m gives a scalar measure of the tendency of the particles to align with 

r— »oo 
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Figure 3: (Color online) Reduced pressure p = PV/NkT in a system of N = 5000 superellipses with 
exponent £ = 7.5 during MD runs with 7 = ±2.5 • 10~ 5 . The predictions of simple versions of SPT 
and FV theory are also shown for comparison. The agreement with FV predictions is not perfect; 
a numerical fit produces a coefficient 2.9 instead of 3 in the numerator of Eq. (|2*|). Particularly 
noticeable are the change in slope around (j) ~ 0.72 and also the transition onto a solid branch 
well-described by free-volume approximation around (j) ~ 0.8. Starting the decompression from an 
ordered tiling in which all rectangles are aligned produces identical pressure to within the accuracy 
available. Systems of N = 1250 and N = 10000 particles, as well as a wide range of particle growth 
rates, were investigated to ensure that there were no strong finite-size or hysteresis effects. In faster 
compressions of an isotropic liquid one gets smaller final densities due to the occurrence of defects 
such as vacancies or grain boundaries. 

a global coordinate system; S 2 is the usual nematic order parameter, and S4 is the tetratic 
order parameter. They can be very easily calculated from an alternative definition 

S m = max (cos[m(# — O )]) , (4) 

#0 
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which can be converted into an eigenvalue problem (in any dimension) for the case m = 2. 20 
When m = 4, we can rewrite it in the same form as m = 2 by replacing 9 with 26. The 
vector n m = (cos 8q, si 11 $o) determines a natural coordinate system for orient at ionally ordered 
phases. It is commonly called the director for nematic phases (m = 2), and we will refer to 
it as a bidirector for tetratic phases (m = 4). 

In two-dimensional liquid-crystalline phases, it is expected that there can be no long- 
range orientational ordering, but rather only quasi-long-range orientational ordering.— Based 
on elasticity theory with a single renormalized Frank's constant K = irK/(8kBT), it is 
predicted^ 2 , that there will be a power-law decay of the correlations at large distances, 
Gm{ r ) ~ r _r? > where 

7] = m 2 /16K. (5) 

This would imply that S m vanishes with increasing system size, 

S m ~ N~^\ (6) 

We note that this prediction is based on literature for the nematic phase. We are not aware 
of any theoretical work explicitly for a tetratic phase. 

The KTHNY theories predict that the isotropic liquid first undergoes a defect-mediated 
second-order transition into an orientationally quasi-ordered but translationally disordered 
state when K = 1 by disclination pair binding. At higher densities there is another second- 
order phase transition into a solid that has long-range orientational order and quasi-long 
range translational order, mediated by dislocation pair binding. The validity of this theory 
is still contested even for hard disks,— and its applicability to systems where there is strong 
coupling between orientational and translational molecular degrees of freedom is question- 
able. Additionally, the basic theory needs to be modified to include three independent elastic 
moduli as opposed to only two in the case of six-fold rotational symmetry. 

The observed change in S4 as an isotropic liquid is slowly compressed is shown in Fig. 
0] for both MD and MC runs. It is clearly seen that tetratic order appears in the system 
around <p ~ 0.7 and increases sharply as the density is increased, approaching perfect order 
(S4 = 1) at close packing. Throughout this run S2 remains close to zero and thus no 
spontaneous nematic ordering is observed. It is important to note that superellipsoids are 
not perfect rectangles and have rounded sides. It is therefore not unexpected that they show 
less of a tendency toward tetratic (right-angle) ordering, and have the isotropic-tetratic (IT) 
transition at slightly higher densities. Additionally, the MD runs show more (correlated) 
variability due to the strong correlations between successive states (snapshots). Therefore, 
we prefer to consider the MC results, other than at very high densities when we have to 
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Figure 4: (Color online) Values of the order metrics S4 [c.f. Eq. @] and [c.f. Eq. ©] for 
snapshot configurations along (de)compression MC and MD runs with N = 5000 particles. Also 
shown is an MC run with N = 1250 for comparison. A transition in S4 is visible around <p ~ 0.7, 
and around </> ~ 0.8 for T^. The values of are much smaller and variable for compression runs 
due to the sensitivity to the exact axes used, however, zooming in reveals a qualitative change in 
Tk around (j) ~ 0.8 even in the compression data. 

resort to MD studies. We have also performed runs decreasing the density of a random 
domino tiling, which has no nematic but has perfect tetratic order, and the resulting £4 is 
also shown in the figure. Only a mild hysteresis is seen, especially for the MC runs, which 
would be indicative of a continuous IT transition, or at least a weakly discontinuous one. 

Figure El shows G^r) for a collection of states in the vicinity of the IT transition, thor- 
oughly equilibriated using MC, on both a log-linear (lower densities) and a log-log (higher 
densities) scale. It is seen that there is a clear change in the long-range behavior of G^(r) as 
the density crosses above C ~ 0.70, from an exponential decay typical of an isotropic liquid, 
to a slower-than exponential decay at higher densities. The decay tails at higher densities 
are rather consistent with a power-law decay, and the fitted exponents rj are shown in Fig. H3 
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Figure 5: (Color online) Left: Log-linear plot of G±(x) for thoroughly equilibriated samples of 
N = 10000 particles, showing the decay of orientational ordering with distance x = r/D. The 
isotropic-tetratic transition occurs between <p = 0.69 and 0.70, when the tail behavior of G±(r) 
changes from exponential (short-ranged) to slower-than-exponential. Right: Log- log plot of G&(x) 
for equilibrated systems of N = 5000 particles, showing power-law decay indicative of quasi- long- 
range tetratic order. The fitted values of the power-law exponent are shown in Fig. H3 

It can be seen that rj crosses the value r] c = 1 predicted by KTHNY theory when »s 0.71, 
which is very consistent with the estimates of the location of the IT transition through the 
other methods above. It is not clear to us why the authors of Ref. loused the value of the ex- 
ponent predicted by KTHNY theory for the bond-bond orientational order in the hard-disk 
system, rj c = 1/4, instead of r] c = 1. The somewhat higher values for S4 for the system with 
iV = 1250 relative to the system with iV = 5000 particles are quantitatively well-explained 
by Eq. (JHJ) using the values of 77 from Fig. El We note that we have never observed a phase 
boundary between a crystallized region and a disordered liquid, which would be indicative 
of a first-order phase transition. 

C. Translational Order 

Measuring translational order is more difficult than orientational ordering. From the 
observations above we are motivated to look for translational order of the kind present in a 
random domino tiling. Looking at the centroids of the dominos themselves does not reveal a 
simple pattern. However, if we split each rectangle into two squares and look at the centroids 
of the 2N squares, translational order will be manifested through the appearance of square- 
lattice periodicity. Such periodicity is most easily quantified by the Fourier transform of the 
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Figure 6: Log-linear plot of I/77, where 77 is the exponent of decay of G±{r) found by fitting the 
G±(x) data in Fig. 0to a power-law curve, G±(x) = Cx~ r '. 

square centroids, i.e., the structure factor 



S(k) 



1 

N 



N 



exp (ik • rj 



(7) 



In a translationally disordered state, 5(/c) is of order one and decays to unity for large k. For 
long-ranged periodic systems, S(k) shows sharp Bragg peaks at the reciprocal lattice vectors, 
while for quasi-long-range order the peaks have power-law wings. It is however difficult to 
exactly determine when true peaks replace the finite humps that exist due to short-range 
translational ordering in the liquid state. 

It would be convenient to have a scalar metric of translational order similar to S4 for 
tetratic order. We use the averaged value of S(k) over the first four Bragg peaks 



1 

2iV 



S(^nn)+£(^n ± ; 
a a 



where a = a/y/4> is the expected spacing of the underlying square lattice, and ny and nj_ 
are two perpendicular unit vectors determining the orientation of the square lattice. In the 
tiling limit Tk = 1, and for a liquid ~ 0. When decompressing a prepared tiling, we 
already know nn = (1, 0) and it is best to use this known value. However, when compressing 
a liquid, we have no way of knowing the final orientation of the lattice and therefore we use 
the bidirector nn = n 4 , as determined during the measurement of 5*4. This method seems 
not to work well because even small fluctuations in the director cause large fluctuations in 
Tk, and additionally, small defects can disrupt periodicity and significantly reduce the value 
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of Tfc below unity. In Fig. 0]we show the values of Tk along with S4. It is seen that for the 
decompression run, Tk starts at unity and decays continuously until it apparently goes to 
zero around ~ 0.8. The compression runs similarly show a visible but noisy and masked 
increase above zero when the density increases above ~ 0.8, but do not reach the same 
level of Tk as for decompression from a perfect crystal. We are therefore led to believe that 
there is a second transition from tetratic liquid to tetratic solid at <fi ~ 0.8. 

In addition to reciprocal space S(k), one can also look at the center-to-center-distance 
distribution function g 2 (r) for the squares (half dominos). However, quantitative analysis 
of gi{r) is made difficult because of oscillations due to exclusion effects and also due to the 
coupling to orientation. Instead of presenting such a one-dimensional pair correlation func- 
tion, we present g 2 (Ax, Ay), which is simply the orientationally-averaged g 2 (Ax, Ay, A6). 
In Figs. [7[|Hland|niwe show a snapshot of a system of iV = 5000 dominos, along with the cor- 
responding g 2 (Ax, Ay) and S(k), for three densities, corresponding to an isotropic liquid, a 
tetratic liquid [i.e., a state with (quasi-) long range tetratic but only short-range translational 
order], and a tetratic solid [i.e., a state with (quasi-) long range tetratic and translational 
order]. For the g 2 (Ax, Ay) plots, we have drawn the expected underlying square lattice at 
that density. Note that g 2 (Ax, Ay) always has two sharp peaks corresponding to the square 
glued to the one under consideration in the dimer (domino). 

D. Solid Phase 

Having confirmed the appearance of a translationally-ordered tetratic phase closely re- 
lated to domino tilings of the plane, we turn to understanding the nature of the the 
thermodynamically-favored tiling. There are two likely possibilities: The tiling shows (trans- 
lational) ordering itself, or the tiling is "random". In the context of a discrete system like 
domino tilings, the concept of a random tiling is mathematically well-defined in terms of 
maximizing entrop y 1 ^ 24 This random tiling has a positive degeneracy entropy 0.58313/cb, 
unlike ordered tilings such as the nematic tiling (in which all dominos are alligned). 

Our compressions of isotropic liquids have invariably led to apparently disordered domino 
tilings upon spontaneous "freezing", albeit with some frozen defects. This suggests that the 
disordered tiling has lower free energy than ordered tilings. However, it is also possible that 
the disorder is simply dynamically trapped when the tetratic liquid freezes. In fact, starting 
a decompression run from an aligned nematic tiling shows that the tiling configuration is 
preserved until melting into a tetratic liquid occurs around <fi ~ 0.8. This is demonstrated in 
Fig. Uni where both S4 and S 2 as well as T k are shown along a decompression run starting 
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Figure 7: (Color online) A snapshot configuration of a system of N = 5000 dominos at (p = 0.7 
(top) with inset with threefold magnification showing local packing structure, along with g<i (Ax, Ay) 
overlayed over the underlying square lattice (bottom left) and S(k) (bottom right), obtained after 
splitting each domino into two squares. It is clear that the system is isotropic from the rotational 
symmetry of S(k). Only short-range order is visible in gz(Ax,Ay), confirming that this is an 
isotropic liquid. 
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Figure 8: (Color online) A system of N = 5000 dominos as in Fig. [7] but at <f> = 0.750, which 
shows a tetratic liquid phase. Four-fold broken symmetry is seen in S(k), but without pronounced 
sharp peaks. The range of ordering in gi (r) has increased, but still appears of much shorter range 
than the size of the system, as seen clearly in the plot of the actual domino configuration. It is 
interesting that g2 (Ax, Ay) is very anisotropic, being much stronger to the side of a square relative 
to its diagonals. No phase boundary characteristic of first-order transitions is visible. 
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Figure 9: (Color online) A system of N = 5000 dominos as in Figs. [7] and |H] but at <f> = 0.825. 
The structure factor shows sharp peaks (maximum value is above 10) on the sites of a (reciprocal) 
square lattice, and g2 (r) shows longer-ranged translational ordering, indicating a solid phase. Visual 
inspection of the configuration confirms that the translational ordering spans the system size and 
shows some vacancies consisting of only a single square (half a particle). 
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Figure 10: Nematic, tetratic, and translational order metrics as a domino tiling in which all rect- 
angles are aligned is slowly decompressed from close-packing. The nematic crystal spontaneously 
realigned to a different orientation of the director from the starting one at around 4> ~ 0.84, causing 
some fluctuations and a drop in which are likely just a finite-size (boundary) effect. 

with both a disordered and an ordered tiling. It is seen that S 2 drops sharply around ~ 0.8 
while 5*4 remains positive until <p « 0.7, clearly demonstrating the thermodynamic stability 
of the tetratic liquid phase in the intermediate density range. Subsequent compression of 
this liquid would lead to a disordered tiling without any trace of the initial nematic ordering. 

We expect that the free-volume contribution to the free energy is minimized for ordered 
tilings at high densities. However, we also expect that solid phase is ergodic in the sense 
that transitions between alternative tiling configurations will occur in long runs of very 
large systems, so that in the thermodynamic limit the space of all tilings will be explored. 
This amounts to a positive contribution to the entropy of the disordered tiling due to its 
degeneracy, and it is this entropy that can thermodynamically stabilize the disordered tiling 



even in the c 



in Rcfs. I5H25I is necessary In particular, including collective Montr Carlo trial moves that 



ose-packed limit. A closer analysis similar to that carried for hard-disk dimers 
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transition between different tiling configurations, as well as relaxation of the dimensions of 
the unit cell (important for smaller solid systems), is important. We believe that, just like 
the hard-disk dimer system, the hard-square dimer system has a thermodynamically stable 
nonperiodic solid phase. Also, as in the hard dumbbell (fused hard-disk dimers) system, 
we expect that for aspect ratios close to, but not exactly, two, the nonperiodic solid will 
be replaced by a nematic (and possibly periodic) phase at the highest densities . 25 > 26 This 
is because reaching the maximal density = 1 seems to require aligning the rectangles. It 
is interesting, however, that at least for rational, and certainly for integer aspect rations 
such as a = 3, there is the possibility of disordered solid phases being stable even in the 
close-packed limit. On physical grounds we expect the phase diagram to vary smoothly with 
aspect ratio, rather than depending sensitively on the exact value of a. 

Accepting for a moment the existence of a nonperiodic solid phase, it remains to verify 
that the compressed systems we obtain in our simulations are indeed like (maximal entropy) 
random tilings of the plane with dimers. This is hard to do rigorously, as it requires compar- 
ing all correlation functions between a random tiling and our compressed systems. Figure 
ITT1 shows a visual comparison of a random tiling of a large square, generated using random 
spanning trees by a program provided to us by the authors of Ref. Q, and a system of 
superellipses compressed to <fi = 0.95 (close to the achievable maximum for our MD program 
for such high superellipse exponents). While the translational ordering in the compressed 
solid is clearly not perfect as it is for the true tiling, visual inspection suggests close similarity 
between the the local tiling patterns of the two systems. In Fig. El we show g2(Ax, Ay) for 
the true tiling, along with the difference in g2 between the true tiling and the compressed 
solid. Here we do not split the rectangles into two squares, i.e., the figure shows the prob- 
ability density of observing a centroid of another rectangle at (Ax, Ay) given a rectangle 
at the origin oriented with the long side along the x axis. It can be seen that there is a 
close match between the random tiling and the compressed solid, at least at the two-body 
correlation level. 

IV. CONCLUSIONS AND FUTURE DIRECTIONS 

The results presented in this paper highlight the unusual properties of the simple hard- 
rectangle system when the aspect ratio is a = 2, hopefully stimulating further research into 
the hard-rectangle system. For square dimers (dominos), in addition to the expected low- 
density isotropic liquid phase, a stable tetratic liquid phase is clearly observed, in which 
there is four-fold orientational ordering but no translational ordering. A tetratic solid phase 
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Figure 11: A comparison between a true random tiling of a square with dominos^ (top), and 
the unit cell of a system of N = 5000 superellipses with exponent £ = 7.5 slowly compressed from 
isotropic liquid to <p = 0.95 (bottom). The compressed system is not a perfect tiling due its lower 
density and frozen defects, as well as the rounding of the superellipses relative to true rectangles. 
Therefore at large scales the two systems look different. However a closer local examination reveals 
similar tiling patterns in the two systems, typical of "random" tilings. 



20 



Figure 12: (Color online) Left: Center-center pair correlation function g2(Ax,Ay) for the perfect 
random tiling in Fig. 1111 This g2 is a collection of ^-functions whose heights can also be calculated 
exactly 27 (the calculation is nontrivial and we have not performed it). We have normalized 52 so 
that the highest peaks have a value of one. Right: The absolute value of the difference between 
$2 (Ax, Ay) for the two systems shown in Fig. Ill[ shown on a coarse-enough scale so that the 
broadening of the peaks due to thermal motion is not visible. The color table used in this figure is 
discrete in order to highlight the symmetry and hide small fluctuations due to finite system size. 
The difference in 52 is almost entirely within the smallest interval of the color table (less than 0.1, 
gray), with only some peaks showing differences up to 0.25. 

closely connected to random domino tilings is found and we conjecture that it is thermody- 
namically stabilized by its positive degeneracy entropy. The transitions between the phases 
are consistent with a KTHNY-like sequence of two continuous transitions. If this is indeed 
the case, then the hard dimer system provides an excellent model for the study of contin- 
uous transitions, with a rather large gap in density between the two presumed transitions 
A(j) « 0.1, unlike the hard-disk system. Random jammed packings of rectangles seem to 
be translationally ordered, similar to the behavior for disks 28 but unlike spheres which can 
jam in disordered configurations.— However, unlike disks, the systems of rectangles show 
orientational disorder, once again illustrating the geometric richness of even the simplest 
hard-particle models. 

Further investigations are needed for the domino system to conclusively determine its 
phase behavior. Improved MC with collective moves that explore multiple tilings, as well as 
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allow for relaxation of the boundary conditions, should be implemented. Additionally, the 
free energies of the different phases should be computed so that the exact locations of the 
phase transitions could be identified. The final goal is to completely characterize the phase 
diagram of the hard rectangle system in the a — <fi plane, as has been done, for example, 
for diskorect angles^ and ellipses.— In addition to nematic and smectic phases, novel liquid 
crystal phases with tetratic order may be discovered. 
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